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ABSTRACT 

Context. The dynamical ejection of stars from star clusters affects the shape of the stellar mass function (MF) in these 
clusters, because the escape probability of a star depends on its mass. This is found in iV-body simulations and has been 
approximated in analytical cluster models by fitting the evolution of the MF. Both approaches are naturally restricted 
to the set of boundary conditions for which the simulations were performed. 

Aims. The objective of this paper is to provide and to apply a simple physical model for the evolution of the MF in 
star clusters for a large range of the parameter space. It should also offer a new perspective on the results from iV-body 
simulations. 

Methods. A simple, physically self-contained model for the evolution of the stellar MF in star clusters is derived from 
the basic principles of two-body encounters and energy considerations. It is independent of the adopted mass loss rate 
or initial mass function (IMF), and contains stellar evolution, stellar remnant retention, dynamical dissolution in a tidal 
field, and mass segregation. 

Results. The MF evolution in star clusters depends on the disruption time, remnant retention fraction, initial-final 
stellar mass relation, and IMF. Low-mass stars are preferentially ejected after t ~ 400 Myr. Before that time, masses 
around 15 — 20% of the maximum stellar mass are lost due to their rapid two-body relaxation with the massive stars 
that still exist at young ages. The degree of low-mass star depletion grows for increasing disruption times, but can be 
quenched when a large fraction of massive remnants is retained. The highly depleted MFs of certain Galactic globular 
clusters are explained by the enhanced low-mass star depletion that occurs for low remnant retention fractions. Unless 
the retention fraction is exceptionally large, dynamical evolution always decreases the mass-to-light ratio. The retention 
of black holes reduces the fraction of the cluster mass in remnants because white dwarfs and neutron stars have masses 
that are efficiently ejected by black holes. 

Conclusions. The modeled evolution of the MF is consistent with iV-body simulations when adopting identical boundary 
conditions. However, it is found that the results from iV-body simulations only hold for their specific boundary conditions 
and should not be generalised to all clusters. It is concluded that the model provides an efficient method do understand 
the evolution of the stellar MF in star clusters under widely varying conditions. 

Key words, stellar dynamics - stars: kinematics - (Galaxy:) globular clusters: general - (Galaxy:) open clusters and 
associations: general - galaxies: star clusters - galaxies: stellar content 



. 1. Introduction a pronounced effect on the inferred globular cluster mass 

y\ . . function (Kruijssen & Portegies Zwart 2009). 

!h The evaporation of star clusters is known to change the . . . . . . . 

d ■ shape of the underlying stellar mass function 1 (Henon The existm ? Parametrised cluster models that mcor- 

1969; Chernoff & Weinberg 1990; Vesperini & Heggie 1997; P orate \ a description of low-mass star depletion are re- 

Takahashi & Portegies Zwart 2000; Portegies Zwart et al. s * ncted h J tb f Physically self-contained models on which 

2001; Baumgardt & Makino 2003). This phenomenon has are h ^J° me studl6S ( LamerS et f L 2006; Kruijssen 

been used to explain the observed MFs in globular clus- ^ La ™ ers 2008) assume an increasing lower stellar mass 

ters (Richer et al. 1991; De Marchi et al. 2007; De Marchi 1 ™J°*? CC °T \^ S . l> u f S - ( , T t 

& Pulone 2007), which are flatter than typical initial mass 2009) fit a changing MF slope to JV-body simulations. In 

c 4-- /Tiifp c i . mcc onmi t i both cases, the models are accurate lor a certain range oi 

functions (IMFs, e.g. balpeter 1955; Kroupa 2001 . In ad- . . , , , , • , , , . , 

dition the effect of a chaneine MF on cluster Dhotome- boundai T conditions, but they do not include a physical 

try has been investigated (Lamers et al. 2006; Kruijssen mod f and are therefore lacking flexibility While A^-body 

& Lamers 2008; Anders et al. 2009). This has been shown ™ latl0n s do mcmde the appropriate physics, they are 

to explain the low mass-to-light ratios of globular clusters very time-consuming As a result, only a limited number of 

(Kruijssen 2008; Kruijssen & Mieske 2009) and to have c usters cai \ be ™lated and the applicability of the sim- 

ulations is thus restricted to the specihc set ol boundary 
conditions for which they have been run. 



* The models presented in this paper will be publicly 

available in electronic form at the CDS via anonymous ftp It would be desirable to obtain a simple physical model 

to http://cdsweb.u-strasbg.fr/ (130.79.128.5) and also at for the evolution of the MF, which would have a short run- 

http://www.astro.uu.nl/~kruijs. time and could be used independently of Af-body simula- 

1 Hereafter, 'mass function' is referred to as 'MF'. tions. Forty years ago, a pioneering first approach to such 
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a model was made by Henon (1969), who considered the 
stellar mass-dependent escape rate of stars from star clus- 
ters. However, the applicability of his model was limited 
due to a number of assumptions that influenced the re- 
sults. First of all, Hcnon (1969) assumed that the clusters 
exist in isolation and neglected the tidal field. As a conse- 
quence, the ejection of a star could only occur by a single, 
close encounter and the repeated effect of two-body relax- 
ation was not included. Secondly, the distribution of stars 
was independent of stellar mass, i.e. mass segregation was 
not included. Both mass segregation and the influence of a 
tidal field are observed in real clusters, and can be expected 
to affect the evolution of the MF. 

The aim of this paper is to derive a physical description 
of the evolution of the stellar MF in star clusters, allevi- 
ating the assumptions that were made by Hcnon (1969). 
This should explain the results found in A-body simula- 
tions and observations, while providing the required flexi- 
bility to explore the properties of star clusters with simple, 
physically self-contained models. The outline of this paper 
is as follows. In Sect. 2, total mass evolution of star clusters 
is discussed. A recipe for the evolution of the MF is derived 
in Sect. 3, covering stellar evolution, the retain of stellar 
remnants, dynamical dissolution and mass segregation. The 
model is compared to A-body simulations in Sect. 4. In 
Sect. 5, the model is applied to assess the evolution of the 
MF for different disruption times and remnant retention 
fractions. The consequences for other cluster properties are 
also considered. This paper is concluded with a discussion 
of the results and their implications. 



2. The mass evolution of star clusters 

The mass of star clusters decreases due to stellar evolution 
and dynamical dissolution. This is expressed mathemati- 
cally as 



for Wo = 7. This timescale implies a mass loss rate due to 
dissolution that can be described with the simple relation 



dM _ (dM\ fdM\ 



(1) 



with M the cluster mass, and the subscripts 'ev' and 'dis' 
denoting stellar evolution and dynamical dissolution. The 
contribution of stellar evolution to the mass loss is derived 
from the decrease of the maximum stellar mass with time 
and depends on the adopted stellar evolution model. 

The dynamical evaporation of star clusters is increas- 
ingly well understood. Over the past years it has become 
clear that clusters lose mass on a disruption timescale tdis 
that is proportional to a combination of the half-mass re- 
laxation time t r h and the crossing time t cr as tdis °c t* h t].~ x 
(e.g. Baumgardt 2001; Baumgardt & Makino 2003; Gieles 
& Baumgardt 2008). It is found that x = 0.75—0.80, de- 
pending on the concentration (c = log (r t /r c )) or King pa- 
rameter (Wo) of the cluster (Baumgardt & Makino 2003). 
This proportionality leads to a disruption timescale that 
scales with the present day mass as (Lamers et al. 2005): 



t dis = t M~<, 



(2) 



with M the cluster mass, to the dissolution timescale pa- 
rameter which sets the rapidity of dissolution and depends 
on the cluster environment, and 7 a constant related to x. 
Lamers et al. (2009) find 7 = 0.62 for W = 5 and 7 = 0.70 



dM 



M M 1 ^ 

tdis ^0 



(3) 



which can be integrated for the mass evolution of the cluster 
due to dynamical dissolution. 

The above formulation of the cluster mass evolution was 
extended to include stellar remnants, photometric cluster 
evolution, and a simple description of the MF in the SPACE 
cluster models (Kruijssen & Lamers 2008). Stellar rem- 
nants were accounted for by assuming initial-final mass re- 
lations (similar to Sect. 3.1 of the present paper), while the 
photometric evolution was computed by integrating stel- 
lar isochrones from the Padova group (Bertelli et al. 1994; 
Girardi et al. 2000). The description of low-mass star deple- 
tion followed the simple model from Lamers et al. (2006) in 
which the minimum stellar mass of the MF increases with 
time. 

The present study provides a new description of the evo- 
lution of the MF which is based on fundamental principles, 
and does not depend on the above prescription for the to- 
tal mass evolution. In addition, the latest Padova models 
(Marigo et al. 2008) are incorporated to calculate the pho- 
tometric cluster evolution. These improvements update the 
SPACE cluster models. 

3. The evolution of the stellar mass function 

To describe the evolution of the MF, the effects of stellar 
evolution, stellar remnant production, and dynamical dis- 
solution need to be included. While the focus of this paper 
lies with the effects of dissolution, a proper treatment of 
stellar evolution is essential. This is described first, before 
presenting a model for cluster dissolution. 2 

3.1. Stellar evolution 

The influence of stellar evolution on the MF is twofold. First 
of all, the maximum stellar mass decreases, because at any 
time during cluster evolution the most massive stars reach 
the end of their lives. Secondly, the stellar remnants that 
are created upon the death of these massive stars constitute 
a part of the MF that can only be lost from the cluster by 
dynamical mechanisms. 

The maximum stellar mass in the cluster as a function of 
its age is taken from the Padova 2008 isochrones (Marigo 
et al. 2008) for mctallicitics in the range Z — 0.0001 — 
0.03. The stellar remnant masses m sr are computed from 
their progenitor stellar mass m using initial-final mass re- 
lations. Following Kruijssen & Lamers (2008), for white 
dwarfs (m < 8 M©) the relation from Kalirai et al. (2008) 
is adopted: 



TOwd = 0.109m + 0.394 M , 



(4) 



which holds for all ages that are covered by the Padova 
isochrones. For neutron stars (8 M© < m < 30 M©) the 
relation from Nomoto et al. (1988) is used: 



m ns = 0.03636(m - 8 M©) + 1.02 M©, 



(5) 



The model presented in this paper is independent of the mass 
loss rate and of the form of the IMF Ni (m) , but for explanatory 
purposes a Kroupa (2001) IMF is adopted later on. 
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while for black holes (m > 30 M Q ) a simple relation is 
assumed that is in acceptable agreement with theoreti- 
cally predicted masses of stellar mass black holes (Fryer 
& Kalogera 2001): 



m bh = 0.06(m - 30 M ) + 8.3 M Q . 



(6) 



With these relations, the remnant MF is computed from 
conservation of numbers as 

N SY (m sr ) = / rct , S r(M)iV(m(TO sr ))-j— , (7) 

with sr = {wd, ns,bh} denoting the appropriate remnant 
type, N SI (m SI ) representing its MF, / re t,sr(M) denoting the 
cluster mass-dependent fraction of these remnants that is 
retained after applying kick velocities, and N(m(m ST )) rep- 
resenting the progenitor MF. 

For a given velocity dispersion of remnants, the reten- 
tion fraction of each remnant type depends on the local 
escape velocity v csc , which is related to the potential ip as 
v csc = \ph\). Stellar remnants are predominantly produced 
in the cluster centre in the case of mass segregation, which 
is reached most rapidly for massive stars (see Sect. 3.2). 
For a Plummer (1911, also see Eq. 9) potential this im- 
plies that upon remnant production v csc = ^/2GM/r , 
with G the gravitational constant and ro the Plummer ra- 
dius. Adopting a Maxwellian distribution of velocities that 
is truncated at v csc , it is straightforward to show that 



fret,sr{x) — A 



XQ 



-x 2 /2 



(8) 



where A is a normalisation constant to account for the trun- 
cation of the velocity distribution and x 2 = 2GM/r a 2 1 ., 
with cr 2 r = (Tq + er 2 ick sr denoting the total velocity disper- 
sion of the produced remnant type, which arises from the 
central velocity dispersion in the cluster a 2 = GM/3ro (e.g. 
Heggie & Hut 2003) and the velocity dispersion of the ex- 
erted kick (Tkick- The normalisation constant then follows as 
A = crfV3 - (-3). 

Typical values of the kick velocity dispersion crkick.sr are 
given in literature. White dwarf kicks have recently been 
proposed to be of order Ckick.wd = 4 km s _1 (Davis et al. 
2008; Fregeau et al. 2009). For neutron stars erkick,ns = 
100 km s _1 is adopted, which is a somewhat conservative 
estimate with respect to theory, but it agrees reasonably 
well with observed neutron star numbers in globular clus- 
ters and represents a compromise between single star and 
binary channels (for estimates of the retention fraction and 
discussions of the 'neutron star retention problem' see Lyne 
& Lorimcr 1994; Drukier 1996; Arzoumanian ct al. 2002; 
Pfahl et al. 2002). Gravitational wave recoils are thought to 
exert black hole kicks of order CTkick.bh = 80 km s _1 (Moody 
& Sigurdsson 2009). This value depends on metallicity, but 
for simplicity I assume a single, typical value here. 

The retention fractions following from Eq. 8 are shown 
as a function of cluster mass per unit Plummer radius in 
Fig. 1. This quantity best reflects the retention fraction 
because x 2 tx M/r$ in Eq. 8. Open clusters (with initial 
masses M\ such that typically M-Jr^ < 3 x 10 M Q pc _1 , 
Larsen 2004) do not retain any neutron stars or black holes, 
while globular clusters {M-Jr ~ 3 x 10 4 — TO 6 M pc" 1 , 
Harris 1996) retain 0.1 — 4% of the neutron stars and 0.3 — 
7% of the black holes. These values are in excellent agree- 
ment with other studies (e.g. Pfahl et al. 2002; Moody & 
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Fig. 1. Retention fraction of stellar remnants as a function of 
cluster mass per unit Plummer radius M/ro, for black holes 
(solid), neutron stars (dashed) and white dwarfs (dotted). 



Sigurdsson 2009), but are still lower than the large observed 
number of neutron stars in a number of globular clusters 
(the aforementioned 'retention problem'). 

3.2. Dissolution and the evolution of the mass function 

Dissolution alters the shape of the stellar MF in star clus- 
ters due to the effects of two-body relaxation and energy 
equipartition. In a pioneering paper, Henon (1969) derived 
the escape rate of stars of different masses from an isolated 
cluster. The cluster was represented by a Plummer (1911) 
gravitational potential: 



V>(r) = Vo 1 + 



-1/2 



(9) 



where r denotes the Plummer radius setting the concen- 
tration of the cluster and tpo = GM/ro represents the cen- 
tral potential, with G the gravitational constant and M the 
cluster mass. It was argued by Henon (1960) that the only 
way for stars to escape such an isolated cluster is by a single, 
close encounter. The corresponding stellar mass-dependent 
escape rate was found to be (Henon 1969): 



dN(m) \E\ 3 / 2 N(m) 
dt ~ GM 9 / 2 



J™ N{m')F (^) m' 2 dm', 



(10) 



with N(m) the MF, m the stellar mass, E the total energy 
of the cluster, and F(p) a function related to the ejection 
probability for a star of mass m in a close encounter with a 
star of mass m' and a corresponding mass ratio fj, = m/m'. 
The expression in Eq. 10 is independent of the adopted 
IMF. The function F will be referred to as the 'Henon func- 
tion' and is shown in Fig. 2. The original expression consists 
of several integrals that have to be solved numerically. In 
Henon (1969), a table is given for the Henon function, but 
it can also be fitted by: 



F{jj) = (0.32 + 0.55/- 35 + 13.26^ 25 ) 1 . 



(11) 



This approaches the power law F{^l) — 0.075398/1 5 / 2 for 
fi > 1, as was derived explicitly by Henon (1969). 
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100.00 



Fig. 2. Henon function F(jj,), which is a measure for the ejection 
probability of a star of mass m in a two-body interaction with 
mass ratio n = m/m' . The dotted line shows the fit from Eq. 11. 



The total mass loss rate corresponding to Eq. 10 con- 
flicts with A^-body simulations (as was already noted by 
Wielen 1971) because only ejections by single, close en- 
counters are included. This restriction implies that the dis- 
ruption timescale tdis is proportional to the crossing time 
(idis, Henon <x t CT ), while A^-body simulations show that it 
scales with a combination of the half-mass relaxation time 
and the crossing time (idis oc ^h^cr 25 ) due to two-body 
relaxation, i.e. the repeated effect of soft encounters (e.g. 
Baumgardt & Makino 2003). Nonetheless, the escape rate 
from Henon (1969) does accurately describe what happens 
if two stars interact and can therefore be used as a start- 
ing point for a more complete description of the evolu- 
tion of the MF. For that purpose, it is convenient to scale 
Eq. 10 to the mass loss rate found in A^-body simulations 
and only use the relative or 'differential' stellar mass de- 
pendence from Henon (1969). This is allowed if the ratio 
idisAdis, Henon = (*rh/icr) ' 75 only depends on global clus- 
ter properties. It is straightforward to show (e.g. Spitzer 
1987; Heggie & Hut 2003) that indeed this is the case as 
irh/icr oc Ay In A with In A the Coulomb logarithm. As 
such, one can write 



dN(m) fdM\ 

— = Ur) dl , *<•»>■ 



(12) 



with (dM/d£)dis the mass loss rate from Eq. 3 (Lamers et al. 
2005) and x(m) the stellar mass-dependent escape rate per 
unit mass loss rate. The quantity x( m ) is completely inde- 
pendent of the prescription for the total mass evolution. In 
order to derive x( m ); I start from Eq. 10 and express x(m) 
as 



X(m) = 



(13) 



N(m) J °° N(m')F(m/m')\(m, m')m' 2 dm' 
J °° m"N(m") J °° N(m')F(m" /m')X{m" , m')m' 2 dm'dm" ' 

where A(m, m') represents a correction factor to account 
for additional physics (see below). The numerator reflects 
the escape rate, while the denominator is proportional to 
the mass loss rate. 



For mathematical simplicity 3 Henon (1969) made the 
following assumptions in the derivation of Eq. 10. 

(1) The cluster exists in isolation and the tidal field is ne- 
glected. Therefore, ejection can only occur by a single, 
close encounter and the repeated effect of soft encoun- 
ters (two-body relaxation) is not accounted for. This 
underestimates the escape rate of massive stars. 

(2) The distribution of stars is independent of stellar mass, 
i.e. mass segregation is not included. Depending on the 
balance between their encounter rate and their proxim- 
ity to the escape energy, this over- or underestimates 
the escape rate of low-mass stars from Henon (1969). 
Considering the results from Baumgardt & Makino 
(2003), the latter seems to be the case. 

The remainder of this section concerns the derivation of 
the factor A(m, m') in Eq. 13 that corrects for the above 
assumptions. 

Let us assume that the distribution of stars over radius 
and velocity space is initially independent of their mass. 
This implies that mass segregation is dynamically created 
and not primordial, which is discussed in Sect. 6. For such 
an initial distribution, the separation from the escape en- 
ergy AE is independent of mass. As the cluster evolves, 
energy equipartition is reached between the stars and the 
radius, velocity and proximity to the escape energy become 
a function of stellar mass. I first consider this effect on the 
escape rate before including the timescale on which two- 
body relaxation occurs for different stellar masses. Please 
note that the formulation of Eq. 13 with A(m, m') appear- 
ing in the numerator and the denominator implies that only 
the proportionality of A(m, m') is important. Its exact value 
is determined by constants that drop out when substituting 
in Eq. 13. 

It is intuitive to express the dependence of the es- 
cape rate on the energy needed for escape as dN(m)/dt oc 
[AE(m)] . The energy that is required for escape AE is 
related to the position and velocity of the star. 4 For the 
potential in Eq. 9 it is given by 



AE(r,v)=iP(r)- — =ip \ I 



-1/2 



V 

T 



(14) 



with r and v the radial position and velocity of the star, 
and u csc = y/2ip(r) its escape velocity. If the cluster is 
in 'perfect' energy equipartition and correspondingly per- 
fect mass segregation, the radius and velocity become a 
monotonous function of stellar mass (Heggie & Hut 2003, 
Ch. 16). Mass segregation is strongest in the cluster centre, 
which for a Plummer (1911) potential can be approximated 
with a harmonic potential ip oc r 2 . For a cluster in a tidal 
field the potential is truncated, and the harmonic approxi- 
mation serves as a crude but reasonable approximation for 



And because this is the only way to obtain an analytical 
solution as in Eq. 10. 

4 The energy difference AE that is discussed here concerns 
the energy that needs to be added to reach the escape energy. 
As such, it differs from the separation from the escape energy in 
Fukushige & Heggie (2000) and Baumgardt (2001), who are con- 
sidering the excess energy of stars and its relation to the escape 
time, resulting in the aformentioned relation idis oc t°j? B ta 2B - 
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the entire cluster (Heggie & Hut 2003, Ch. 16). Energy 
equipartition yields 



2/ \ / \2 \ m ) 

v (to) = (v) , 



(15) 



with 



oc tpo the mean speed of all stars squared and 



(to) the mean stellar mass. For the harmonic potential, this 
translates to a similar relation for the radial position: 



r(m) = ro 



to 



(16) 



where ro represents the typical radius of the system, in 
this case the Plummer radius. This relation assumes that 
there is no particular stellar mass which dominates the mass 
spectrum. The decrease of radial position with stellar mass 
implied by Eq. 16 is a direct consequence of the energy loss 
endured by massive stars 5 as the system evolves towards en- 
ergy equipartition. Substituting Eqs. 15 and 16 into Eq. 14 
and dividing out the proportionality (v) 2 oc ipo gives an 
expression for AE(m): 



AE(m) = 1 + 



(m) 



-1/2 



ci- 



(to) 



(17) 



with ci = (v) 2 /2tp denoting the ratio of the mean speed 
squared to the central escape velocity squared. This con- 
stant mainly depends on the degree of mass segregation. 
Consequently, it will depend on the IMF. By comparing 
the models to the iV-body simulations with a mass spec- 
trum by Baumgardt & Makino (2003) the value is con- 
strained to ci = 0.020 for a Kroupa IMF, using King (1966) 
potentials with King parameter Wo = 5 — 7 (see Sect. 4). 
For reference, an unevolved Plummer (1911) potential has 
(v) 2 /2ip = 3tt/64 = 0.147. 

By comparing the models to iV-body simulations (pro- 
vided by M. Gieles, private communication) with different 
IMF power law slopes a and a ratio between the maxi- 
mum and minimum mass of 10, the approximate relation 
logci w a — 3.76 is found 6 for a MF n s oc m~ Q . Fitting 
the Kroupa IMF with a single power law in the mass range 
0.08—15 M (as used by Baumgardt & Makino 2003) yields 
a = 2.06, resulting in c\ = 0.020 as mentioned earlier. 7 
The comparison with TV-body simulations also showed that 
a single value of c\ suffices to determine the MF evolution, 
even though it does not remain constant over the full cluster 
lifetime. 

Because dN(m)/dt oc [A.E(m)] _1 , Eq. 17 indicates that 
the escape rate of low-mass stars is increased if a cluster is 
in complete energy equipartition. However, the timescale on 
which two-body relaxation occurs between different stellar 
masses has not yet been considered. For a cluster starting 
with a stellar mass-independent distribution of radial po- 
sitions and velocities, the equipartition timescale t e scales 
as 



t e (m, to') oc to 1 



Til 



1-1 



(18) 



5 And the energy gain experienced by low-mass stars. 

6 This prescription for ci implies that the condition for the 
stars in the cluster to be physically bound AE(m) > is satis- 
fied for all a < 3.63. 

7 Nonetheless, the relation for ci should be expected to exhibit 
some variation for different mass ranges. 



for equipartition between stars of masses m and to,' (Heggie 
& Hut 2003). This is a modified version of the relaxation 
timescale, which shows a very similar proportionality (t T oc 
to -2 ). It illustrates that two-body relaxation occurs on a 
shorter timescale for massive stars than for low-mass stars, 
increasing their escape rate dN(m)/dt oct^ 1 . 

The correction factor for the escape rate A(m, to') that 
appears in the integrals of Eq. 13 now follows from Eqs. 17 
and 18 as 



A(to,to') = t,T 1 (m,m')[A.E(m)]- 1 

-1/2 



= mm 



1 + 



(to) 



TO 



ci- 



(to) 



TO 



(19) 



It was mentioned before that the proportionalities of 
AE(m) and t c (to, to') rather than their exact values suffice 
for the computation of A(to, to') due to the renormalisation 
of the total mass loss rate that appears in Eq. 13: only the 
stellar mass-dependence is important. 

The influence of the tidal field is now included in two 
ways. First of all, the ejection of stars no longer occurs 
by a single, close encounter but arises due to two-body re- 
laxation on the equipartition timescale, representing the re- 
peated effect of soft encounters. Secondly, the above deriva- 
tion of the separation from the escape energy assumes a 
potential which approximates tidally limited clusters. As a 
result, the escape rate of massive stars is increased with re- 
spect to clusters in the model of Henon (1969), which was 
derived for an isolated cluster. On the other hand, the ef- 
fect of mass segregation is included by introducing a stellar 
mass-dependence for the energy needed by stars to reach 
the escape velocity. Low-mass stars are closer to the tidal 
radius than massive stars, leading to a lower energy that is 
needed for escape and an increased escape rate. It depends 
on the shape of the MF which mechanism dominates. 

The evolution of the MF of various cluster components 
is obtained from Eqs. 12, 13 and 19 by writing 



dlogA^ C omp(TO) dlog N(m) 



dt 



dt 



(20) 



where the MFs of stars, white dwarfs, neutron stars and 
black holes are represented by Af comp (m), with comp = 
{s, wd, ns, bh}. The overall cluster evolution is computed by 
combining the results of this section with the prescription 
for stellar evolution from Sect. 3.1. 

If stellar evolution is included, the resulting mass loss 
causes an expansion of the cluster, during which stars are 
lost independently of their masses. This delays the onset 
of mass segregation and the stellar mass-dependent mass 
loss that is described above. The moment of transition to 
stellar mass-dependent mass loss can be characterised by 
a certain fraction of the initial cluster mass that has been 
lost by dissolution /di ss = Mdiss/M- It is assumed that the 
fraction / sm( j of the mass loss for which the ejection rate 
depends on the stellar mass grows exponentially 8 between 
and 1 as 



./; 



smd 



(J I g/diss//diss,scg ^ 



(21) 



8 This form assumes that the increase of the fraction of the 
mass loss that is stellar mass-dependent scales with the total dy- 
namical mass loss, which is a compromise between a step func- 
tion and a linear increase. 
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Fig. 3. MF slope change Aa in the range m = 0.1 — 0.5 Mq 
versus the remaining mass fraction for a Kroupa IMF (solid), 
Salpeter IMF (dotted), and a power law IMF with a = 1.35 
(dashed). In all cases, the IMF mass range is m = 0.1 — 1 M . 
The displayed relation is valid if stellar evolution is excluded. 

where the subscript 'smd' denotes 'stellar mass-dependent', 
/diss, SC g = Afdiss,seg/M is the fraction of the initial mass 
that has been lost by dissolution at which mass segrega- 
tion is reached, and C — (e — 1) _1 is a constant to nor- 
malise / smd = 1 at the reference value / diss = /di ss , sog - For 
/diss > /diss,sc g , per definition / smd = 1, indicating that 
all mass loss is stellar mass-dependent. The timescale i sog 
on which mass segregation is reached and the transition to 
stellar mass-dependent mass loss is completed is propor- 
tional to the initial half-mass relaxation time (t scg oc f r h,i)- 
It has been shown in several studies that for Roche lobe- 
filling clusters the disruption timescale idis oc t^ft^ 25 
(Vesperini & Heggie 1997; Baumgardt & Makino 2003; 
Gieles & Baumgardt 2008), implying that £ scg Adis oc idis 3 - 
The expression for idis in Eq. 2 then leads to i seg /idis oc 
£q- 33 M°' 337 . Assuming that the cluster mass evolution is 
close to linear, the first-order relation /di ss ,sc g oc i sog /idis is 
obtained, implying 

\ 0.33 / ^ 0.21 

[wi^j ' (22) 

for a King parameter of W — 5, with the dissolution 
timescale at the solar galactocentric radius t$ = 21.3 Myr. 
For a King parameter of Wo — 7, the exponent of the initial 
cluster mass M[ becomes 0.23 and tf = 10.7 Myr (Kruijssen 
& Mieske 2009). In this relation, c 2 represents a constant 
that is fixed by comparing the model to the results of TV- 
body simulations from Baumgardt & Makino (2003), giving 
c 2 = 0.25 for W = 5 and c 2 = 0.15 for W = 7 (see Sect. 4). 
The variation with King parameter arises because two-body 
relaxation is faster for more concentrated clusters. If stellar 
evolution were neglected, at all ages c 2 = and / sm d = 1. 

The modeled MF slope change Aa in the mass range 
m = 0.1 — 0.5 M© is shown in Fig. 3 for different IMFs 
covering m = 0.1 — 1 M©. Evidently, Aa is a function of 
the remaining mass fraction and is insensitive to the slope 
of the IMF, as long as that the ratio between the maximum 
and minimum mass is kept fixed and stellar evolution is 
excluded. This is an interesting observation in view of the 



Fig. 4. Relative ejection rate x( m )/N(m) as a function of 
stellar mass, shown for a Kroupa MF with different maximum 
masses. The end point of each curve (dot) marks its maximum 
mass. The quantity x{m)/N(m) = (dlogiV(m)/di)/(dM/di) 
represents the ejection rate per unit mass loss rate normalised 
to the number of stars at each mass (also see Eq. 13). 



MF evolution of globular clusters, in which m w 0.1 — 1 M© 
and stellar evolution only plays a minor role. Figure 3 shows 
that the slope of the MF in globular clusters could be a 
possible indicator for the mass fraction that has been lost 
due to dissolution, provided that the IMF does not vary 
and the remnant retention fractions were not substantially 
dissimilar during the early evolution of different globular 
clusters (see Sect. 5.2 and Fig. 19). 

For the particular example of a Kroupa MF that is trun- 
cated at different maximum masses m mM , the relative es- 
cape rate per unit mass loss rate x(m)/iV(m) (see Eqs. 12 
and 13) is shown in Fig. 4. This quantity is proportional to 
dlog AT(m)/di and reflects the probability that a star of a 
certain mass is ejected. Figure 4 illustrates that the mass of 
the highest relative ejection rate is related to the maximum 
mass of the MF. The peak occurs at intermediate masses if 
the MF is truncated at a high mass. This implies that there 
is a typical mass where the stars are not too far from the 
escape energy and have an equipartition timescale with the 
massive stars that is short enough to eject them efficiently. 
This 'sweet spot' depends on the maximum mass of the 
MF. If the MF is truncated at an intermediate mass, the 
combination of quick two-body relaxation and proximity to 
the escape energy favours the ejection rate of stars at the 
lowest masses. 

The maximum stellar mass at which the transition from 
'sweet spot'-depletion to low-mass star depletion happens, 
is determined by the proximity of the low-mass stars to 
the escape energy. In Fig. 5, the mass of the peak relative 
ejection rate is shown as a function of the maximum stellar 
mass. At low truncation masses, the peak occurs at the 
minimum mass, indicating strong low-mass star depletion. 
Around m max ~ 3 M©, the relative ejection rate at m poa k ~ 
0.4 M© becomes larger than its value at the lowest masses, 
which causes a jump in Fig. 5. For even higher values of 
w max , the peak relative ejection rate typically occurs at 
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Fig. 5. Mass of the highest relative ejection rate m pea k as a 
function of the maximum stellar mass of the MF ra max (solid 
line). The dashed line represents the relation m pcak = 0.2m max , 
while the dotted line describes an eyeball fit for masses m max > 
3 M and includes an exponential truncation at the low-mass 
end (see Eq. 23). 

15 — 20% of the maximum mass, approximately following 
the relation 

m pcak - 0.2m max c- 2M °/ m — . (23) 

Even though its quantitative properties only hold for a 
Kroupa MF, the variation of the relative escape rate with 
the maximum mass of the MF has several implications 
for star cluster evolution. The change of m max in Figs. 4 
and 5 can be interpreted as an example of what happens 
when stellar evolution removes the most massive stars in 
the cluster, provided that the remnants are all ejected by 
their kick velocities. If dynamical evolution does not affect 
the shape of the MF too much before 7n max (i) <~ 3 M Q , or 
t ~ 400 Myr, the subsequent evolution of the MF will be 
dominated by low-mass star depletion. If substantial dis- 
solution occurs earlier on, it is dominated by the 'sweet 
spot' depletion of intermediate masses. Only the retention 
of massive stellar remnants will make the evolution of the 
MF deviate from these basic estimates, because remnant 
retention can provide a fixed maximum (remnant) mass of 
the MF. This is treated in more detail in Sect. 5. 

4. Comparison to iV-body simulations 

The model described in Sect. 3 can be easily verified 
by running it for the exact same boundary conditions as 
the iV-body simulations 9 by Baumgardt & Makino (2003) 
and comparing the results. They conducted simulations of 
Roche lobe- filling clusters between 8k and 128k particles, 
which were evolved in the Galactic tidal field at galactocen- 
tric radii in the range 2.833 — 15 kpc. The boundary condi- 
tions for the A^-body runs of Baumgardt & Makino (2003) 
differ from those described in Sect. 3 by neglecting kick ve- 
locities and defining the Kroupa stellar IMF between 0.1 
and 15 M Q , thereby excluding black holes. For this partic- 
ular comparison, the same IMF, stellar evolution prcscrip- 

9 These were performed using NB0DY4 (Aarseth 1999). 



tion, and initial-final mass relation for stellar remnants are 
used in the model that is presented in this paper. 

In Fig. 6, the modeled evolution of the (luminous) stellar 
MF is compared to the iV-body runs with King parameter 
Wo = 5 for a range of cluster masses and total disrup- 
tion times. As time progresses, the maximum stellar mass 
decreases due to stellar evolution and the MF is lowered 
due to the dynamical dissolution of the star cluster. The 
slope of the MF changes due to the preferential ejection of 
low-mass stars, which have energies closer to their escape 
energies, even to the extent that it dominates over their rel- 
atively slow two-body relaxation. For both the models and 
the A^-body simulations, the MF develops a slight bend at 
m <~ 0.3 M Q when approaching total disruption. The bend 
arises as an optimum between on the one hand high energies 
but slow relaxation for the lowest stellar masses, and on the 
other hand quick relaxation but low energies for the highest 
stellar masses (see the discussion at the end of Sect. 3). 

In all cases, the resemblance of the models and the N- 
body simulations is striking. The models reproduce all key 
aspects of the A-body runs, such as the amount of low- 
mass star depiction, the changing slope at m ~ 0.3 M Q 
for clusters close to dissolution, the survival of the Kroupa 
bend at m — 0.5 M Q , and the dependence of the low-mass 
depletion on the total lifetime of the cluster (compare the 
three 32k runs). The only difference occurs at the high-mass 
end of the MF, where the maximum stellar masses do not 
match at young ages. This is due to a minor dissimilarity 
of the total mass evolution (also see Lamers et al. 2005; 
Kruijssen & Lamers 2008). Because the maximum stellar 
mass only depends on the age of the cluster, this causes a 
difference in maximum stellar mass when showing the MFs 
at fixed remaining cluster mass fractions. The contrast is 
clearest at young ages, since there the maximum stellar 
mass most rapidly decreases. 

In the description of the model in Sect. 3, two con- 
stants have been determined from the A^-body simulations 
by Baumgardt & Makino (2003). These constants are the 
ratio of the mean speed squared to the central escape ve- 
locity squared (ci, see Eq. 17) and the proportionality con- 
stant for the relation marking the transition to stellar mass- 
dependent mass loss (c2, see Eq. 22). As mentioned in 
Sect. 3, for a Kroupa IMF and King parameter Wo = 5 
one obtains c\ = 0.020 and c 2 = 0.25. To illustrate the 
robustness of the models, in Fig. 7 they are compared to 
a 64k A^-body run with Wq = 7. For such a cluster with 
a higher concentration, the early mass segregation implies 
c 2 = 0.15. Again, the model and the simulation are in ex- 
cellent agreement. 

The dependence of the MF evolution on both constants 
is considered in Fig. 8. For ci, the dependence of the evo- 
lution of the MF on its value is shown in the upper panel 
of Fig. 8, while for c 2 it is shown in the bottom panel of 
Fig. 8. Both panels show the evolution of the MF for the 
64k cluster in Fig. 6 for different values of c\ and c 2 . 

The ratio of the mean speed squared to the central es- 
cape velocity squared c\ affects the ejection probability of 
the stars with the lowest masses. Because these stars are 
closest to their escape energies in a mass-segregated cluster, 
they are most strongly influenced by the value of c\. For 
higher c\ , the MF gets more depleted in low- mass stars due 
to their closer proximity to the escape energy, while for 
lower ci more low-mass stars are retained as the balance 
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Fig. 6. Comparison of the evolution of the stellar MF from the models (dashed) to the iV-body runs from Baumgardt & Makino 
(2003, solid) for the exact same boundary conditions. The initial number of particles and the galactocentric radius are indicated 
in the bottom-left corner of each panel. From top to bottom, the subsequent MFs in each panel are shown for the times at which 
the remaining cluster mass fraction equals M/M = {1, 0.75, 0.6, 0.5, 0.3, 0.2, 0.15, 0.1}. 



between close proximity to the escape energy and slow re- 
laxation shifts to the latter. 



The proportionality constant for the transition to stellar 
mass-dependent dissolution c-i in Eq. 22 affects the MF as a 
whole. For lower C2, the transition occurs earlier and more 
low-mass stars are lost, while for higher the onset of the 
depletion is delayed and the slope of the MF remains closer 
to its initial value. If one were to assume a constant /diss.scg, 
which is contrary to the adopted relation with cluster mass 
in Eq. 22, this would therefore yield a stellar MF in massive 
clusters that is underpopulated in low-mass stars, and a MF 
in low-mass clusters that is overabundant in low-mass stars. 



5. Star cluster evolution 

In this section, the described model is applied to compute 
the evolution of clusters for a variety of boundary condi- 
tions. The stellar content as well as integrated photometry 
are addressed, using the boundary conditions from Sect. 3 
instead of those that were adopted to compare the model 
to iV-body simulations in Sect. 4. The most important dif- 
ferences are the mass range of the IMF, the inclusion of 
remnant kick velocities, and the initial-final mass relation. 

The model that will be referred to as the 'standard 
model' uses a metallicity Z — 0.004 (which is typical of 
globular clusters), a King parameter 10 of Wo — 7 (cor- 
responding to 7 = 0.7 in Eq. 2), a dissolution timescale 
parameter to ~ 1 Myr, and a Kroupa IMF between 



For Wo = 5, or 7 = 0.62, the results vary only marginally. 
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Fig. 7. Comparison of the evolution of the stellar MF from 
the models (dashed) to the iV-body run from Baumgardt & 
Makino (2003, solid) with Wo = 7 for the exact same boundary 
conditions. From top to bottom, the subsequent MFs are shown 
for the times at which the remaining cluster mass fraction equals 
M/Mi = {1, 0.75, 0.6, 0.5, 0.3, 0.2, 0.15, 0.1}. 



m = 0.08 Mq and the maximum stellar mass given by 
the Padova isochrones at logt = 6.6, which is typically 
m <~ 70 M©. For the computation of the retained remnant 
fraction (see Eq. 8), the Plummer radius ro is related to the 
half-mass radius as = 1.3ro. The half-mass radius is 
assumed to remain constant during the cluster lifetime (e.g. 
Aarseth & Heggie 1998). For the relation between and 
initial cluster mass Mi the expression from Larson (2004) 
is adopted: 



rh = 3.75 pc 



M 



10 4 M 







o.i 



(24) 



The models that are used in this section are computed from 
10 7 yr to 1. 65 x 10 10 yr (the maximum age of the Padova 
isochrones) for initial masses between 10 2 M Q and 10 7 M©, 
spaced by 0.25 dex intervals. 

5.1. The influence of the disruption time 

The disruption time of a cluster affects the evolution of 
the MF and of the integrated photometric properties. To 
assess the influence of the disruption time on cluster evolu- 
tion, clusters with low and high remnant retention fractions 
should be treated separately, because the presence of mas- 
sive remnants also has a pronounced effect on the results 
(see Sect. 5.2). As shown in Fig. 1, for a given kick veloc- 
ity dispersion the remnant retention fraction is set by the 
cluster mass. This means that the division between low and 
high remnant retention fractions can be made by making a 
cut in initial cluster mass. 

In Fig. 9, the impact of the disruption time on the evo- 
lution of the MF is shown for a cluster with initial mass 
log (Mi/M©) = 4.5, representing the evolution for low rem- 
nant retention fractions. 11 The range of the dissolution 
timcscale parameter to and resulting total disruption times 



11 High remnant retention fractions will be treated in the dis- 
cussion of the influence of the retention fraction in Sect. 5.2. 




m/M;; 



Fig. 8. Influence of the constants ci and C2 on the evo- 
lution of the stellar MF. From top to bottom, the sub- 
sequent MFs in each panel are shown for the times at 
which the remaining cluster mass fraction equals M/Mi = 
{1,0.75,0.6,0.5,0.3,0.2,0.15,0.1}. Top panel: the values ci = 
{0.010,0.020,0.030} are represented by dashed, solid and 
dotted lines, respectively. Bottom panel: the values C2 = 
{0, 0.25, 0.40} are represented by dashed, solid and dotted lines, 
respectively. For both ci and C2, the second (boldfaced) value 
is the one obtained from the comparison to the TV-body simula- 
tions with Wo — 5 in Fig. 6. 



that are considered in Fig. 9 cover two orders of magnitude. 
As the total lifetime increases, the depletion of the low-mass 
stellar MF close to total disruption becomes more promi- 
nent. Conversely, the MF of short-lived clusters is depleted 
around m ~ 1 M©. As introduced in the last paragraphs 
of Sect. 3, this difference is caused by the fixed timescale 
on which stellar evolution decreases the maximum stellar 
mass, implying that the masses of the most massive stars 
are larger in quickly dissolving clusters than in slowly dis- 
solving ones. Because in short-lived clusters the massive 
stars are still present when the bulk of the dissolution oc- 
curs, their rapid two-body relaxation with intermediate- 
mass stars dominates over the relatively close proximity to 
the escape energy of low-mass stars, yielding a depletion 
at intermediate masses. In long-lived clusters, this cannot 
occur because the very massive stars have disappeared be- 
fore the mass loss by dissolution becomes important, thus 
resulting in the depletion of the very low-mass end of the 
MF. As a rule of thumb, for t < 400 Myr (which is the life- 
time of a 3 M© star) the depletion typically occurs around 
15 — 20% of the mass of the most massive star (see Sect. 3). 
In terms of the total disruption time, the transition from 
intermediate-mass star depletion to low-mass star depletion 
1 Gyr. 



occurs around ij°* al 
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Fig. 9. Influence of the disruption time on the evolution 
of the stellar MF for a cluster with a low remnant reten- 
tion fraction (log (M[/Mq) — 4.5). From top to bottom, the 
subsequent MFs in each panel are shown for the times at 
which the remaining cluster mass fraction equals M/Mi — 
{1, 0.75, 0.6, 0.5, 0.3, 0.2, 0.15, 0.1}. 



A quantifiable way to look at the evolution of the stel- 
lar MF in star clusters is to consider the slope of the 
MF n s oc mT a in certain mass intervals (Richer et al. 
1991; De Marchi et al. 2007; De Marchi & Pulone 2007; 
Vesperini et al. 2009). For the commonly used mass inter- 
vals 0.1 < m/M < 0.5 (ai) and 0.3 < m/M < 0.8 (a 2 ), 
Fig. 10 shows the evolution of the slope a for the same 
clusters as before. Like Fig. 9, this illustrates that for short 
disruption times the slope steepens as the cluster dissolves, 
while for long disruption times the slope flattens with time. 
The presented models and other model runs indicate that 
a\ increases with time for ^°* al < 1 Gyr and decreases 
for ijjis > 2 Gyr. For total disruption times in between 
these values, the slope first increases and then decreases. 
The slope in the second mass interval a 2 shows the same 
behaviour. It increases for ^°* al < 0.5 Gyr and decreases 
for t*°* al > 1 Gyr. 



a 1 7 




0.6 0.4 

M/M, 

Fig. 10. Influence of the disruption time on the stellar MF slope 
a in the range 0.1 < m/M Q < 0.5 (solid) and 0.3 < m/M Q < 0.8 
(dashed) for a cluster with a low remnant retention fraction 
(log (Mi/Mq) = 4.5). Shown is a versus the remaining cluster 
mass fraction. From top to bottom, for each mass range the 
lines represent to — {0.1, 1, 10} Myr, corresponding to t^l" 1 — 
{0.16, 1.42, 12.26} Gyr. 
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Fig. 11. Influence of the disruption time on the M/Ly ra- 
tio evolution for a cluster with a low remnant retention frac- 
tion (log (Mi/M ) = 4.5). Shown is the relative M/L v ratio 
decrease with respect to the value expected for stellar evolu- 
tion (M/Ly) st0 v versus the remaining cluster mass fraction. The 
solid, dashed and dotted lines represent to = {0.1,1,10} Myr, 
respectively, corresponding to t^ 1 = {0.16, 1.42, 12.26} Gyr. 



The mass-to-light (M/L) ratio evolution of star clusters 
is affected by the evolution of the MF due to the large 
variations in M/L ratio between stars of different masses. 
Massive stars have lower M/L ratios than low- mass stars, 
implying that a cluster with a MF that is depleted in low- 
mass stars will have a reduced M/L ratio (Baumgardt & 
Makino 2003; Kruijssen 2008; Kruijssen & Lamers 2008). 
As such, one would also expect a correlation between the 
slope of the MF and M/L ratio. 

In Fig. 11, the evolution of the ratio of the V-band 
M/Ly to the mass-to-light ratio due to stellar evolution 
(M/Ly ) stC v is shown for the same clusters as in Figs. 9 
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Fig. 12. Influence of the disruption time on the combined 
evolution of the MF slope a and the M/Ly ratio for a clus- 
ter with a low remnant retention fraction (log (M;/Mq) = 
4.5). Shown is a versus the relative M/Ly ratio decrease due 
to dynamical evolution. All clusters start at the vertical line 
(M / Ly) / (M / Ly) stev = 1- Solid lines denote the slope in the 
mass range 0.1 < m/M < 0.5, dashed lines designate the 
mass range 0.3 < m/M© < 0.8, and dotted lines represent 
the mass range 0.3m max (t) < iti/Mq < 0.8m max (i), with 
from top to bottom to — {0.1, 1, 10} Myr, corresponding to 
4°s al = {0-16, 1.42, 12.26} Gyr. 

and 10. This quantity reflects the relative M/Ly ratio 
change induced by dynamical evolution with respect to evo- 
lutionary fading only. If the ejection rate would be indepen- 
dent of stellar mass, the evolution would follow a horizon- 
tal line at (M/Ly)/(M/Ly) stev = 1. However, when ac- 
counting for dynamical evolution, the M/L ratio is always 
smaller than that for stellar evolution only. Somewhat sur- 
prisingly, this is also the case for clusters for which the 
slope of the MF increases (see Fig. 10). This is explained 
by looking at the evolution of the entire MF in Fig. 9. Even 
though the slope at low masses increases for short disrup- 
tion times due to the ejection of intermediate-mass stars, 
the most massive stars that dominate the cluster light are 
still retained. Because stars of intermediate masses are lost 
instead, the M/L ratio decreases. 

Because the slope of the stellar MF either increases or 
decreases at masses m < 1 M Q , the decrease of the M/L 
ratio implies a large range of MF slopes that can occur 
at low M/L ratios. This is shown in Fig. 12, where the 
relation between a and the M/L ratio drop is presented. 
The slope of the stellar MF in a certain mass range does 
not necessarily reflect the M/L ratio of the entire cluster. 
Considering the aforementioned rule of thumb stating that 
for total disruption times < 1 Gyr the depletion of the 
MF occurs around 15 — 20% of the mass of the most massive 
star m max (t), it is useful to define the slope in a mass range 
that is related to m max (t). In Fig. 12, the relation between 
slope and M/L ratio is also shown for the slope in the stellar 
mass range 0.3m max (i) < m/M Q < 0.8m max (i). In such 
a relative mass range, the slope follows a much narrower 
relation with M/L ratio. The range between 30% and 80% 
of m max (t) was chosen to maximise this effect. 

For the slopes in the fixed stellar mass ranges (ai 
and ct2, see above), the relation with the M/L ratio be- 
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Fig. 13. Influence of the disruption time on the V — I colour for 
a cluster with a low remnant retention fraction (log (Mi/M©) = 
4.5). Shown is the colour offset due to dynamical evolution 
A(V — I) versus the remaining mass fraction. The solid, dashed 
and dotted lines represent to = {0.1,1,10} Myr, respectively, 
corresponding to 4<?* al = {0.16, 1.42, 12.26} Gyr. 

comes better defined for long-lived clusters. It is shown in 
Figs. 10 — 12 that both the slope and the M/L ratio de- 
crease for clusters with long disruption times, indicating 
that both quantities are more clearly related for globular 
cluster-like lifetimes. 

The colour of star clusters is also influenced by the evo- 
lution of the MF, due to the colour differences between 
stars of different masses. The V — I magnitude difference 
A(V — I) with respect to the V — I value that a cluster 
would have if dynamical evolution were neglected is shown 
in Fig. 13. As the clusters dissolve, their colours become 
redder due to the ejection of main sequence stars. The mag- 
nitude difference in V — I exceeds A(V — I) = 0.1 mag for 
total disruption times < 1.5 Gyr. In redder passbands (e.g. 
the V — K colour), the difference grows to several tenths 
of magnitudes. For longer total disruption times only stars 
of the lowest masses are ejected (see Fig. 9), which hardly 
contribute to the cluster light and colour, implying that the 
colours are only marginally affected. 

5.2. The influence of the remnant retention fraction 

The formation of stellar remnants introduces massive bod- 
ies in the MF that do not end their lives due to stellar evo- 
lution like massive stars do. Depending on their kick veloc- 
ities, stellar remnants can be retained in (massive) clusters. 
If they are retained, they keep affecting the evolution of the 
stellar MF until the cluster is disrupted. Especially black 
holes can have a pronounced effect on cluster evolution. 

The remnant retention fraction arises from the cluster 
mass, radius and the kick velocity dispersion (see Eq. 8). In 
this section, the mass-radius relation from Eq. 24 is used. 
Although the results will differ for other relations, it has 
been verified that for commonly used alternatives, 12 the 
change is only marginal and does not affect the nature of 
the conclusions. To separate the effect of remnant reten- 
tion from that of the disruption time, a fixed initial cluster 

1 Such as a constant radius or density. 
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Fig. 14. Influence of the black hole kick velocity dispersion 
and disruption time on the evolution of the stellar MF for 
an initial cluster mass Mi = 10 Mq. From top to bottom, 
the subsequent MFs in each panel are shown for the times 
at which the remaining cluster mass fraction equals M/Mi = 
{1, 0.75, 0.6, 0.5, 0.3, 0.2, 0.15, 0.1}. Solid lines denote to = 1 Myr 
0d°s al = 15-13 Gyr), while dotted lines represent to =0.1 Myr 
= 1.66 Gyr). 



mass of 10 6 M Q is assumed while independently varying the 
velocity dispersion of the remnant kick velocities and the 
disruption time. The corresponding evolution of the stel- 
lar MF is shown in Fig. 14, for the standard model (see 
the beginning of this section) with black hole kick velocity 
dispersions <Tkick,bh = {40,80,200} km s _1 , equivalent to 
/ret,bh = {0.219, 0.041, 0.003} for a 10 6 M Q cluster, and for 
dissolution timescale parameters t = {0.1, 1} Myr, which 
for a 10 6 M cluster implies $?* al = {1.66,15.13} Gyr. 
Assuming an age of 12 Gyr, the present-day mass in the 
case of to = 1 Myr is about M ~ 6 x 10 4 M Q , comparable 
to globular clusters. The remaining fraction of the initial 
mass is M/Mi ~ 0.06. 

If the velocity dispersion of black hole kicks is low and 
a relatively large fraction of black holes is retained, then 



3 



o 1 







1 1 1 1 1 1 


i | i i i | i i 


i ' ' ' 4 






] 


— 




/: 




^"^r^ 


/ ; 


- 




X ~- 

\\ " 






\\ - 






\\ ; 






\ V 










1.0 0.8 0.6 0.4 0.2 0.0 

M/M, 

Fig. 15. Influence of the black hole retention fraction on the 
stellar MF slope a in the range 0.1 < m/M© < 0.5 (solid) 
and 0.3 < m/M Q < 0.8 (dashed) for an initial cluster mass 
Mi = 10 6 Mq. Shown is a versus the remaining cluster mass 
fraction. From top to bottom, for each mass range the lines 
represent <7 kickib h = {40,80,200} km s _1 , corresponding to 
Aet.bh = {0.219,0.041,0.003} for a 10 6 M cluster. 



the ejection rate of massive stars is increased with respect 
to high kick velocity dispersions. This arises due to the 
quick two-body relaxation between the massive stars and 
the black holes, which will have masses larger than the most 
massive stars after a few Myr of stellar evolution. As a re- 
sult, the ejection rate of low-mass stars is largest in clusters 
containing only few black holes. This happens for clusters 
with either long or short disruption times, but the effect is 
largest for long-lived clusters (the solid lines in Fig. 14). In 
these clusters the maximum stellar mass is more strongly 
decreased by stellar evolution than in short-lived clusters, 
implying that the black hole masses are larger compared to 
the most massive stars in these clusters. For long disrup- 
tion times, the presence of massive remnants therefore has a 
more pronounced effect on the ejection rate of massive stars 
than for short disruption times. If these long-lived clusters 
retain a sufficiently large fraction of the stellar remnants, 
their stellar MF may even become depleted in massive stars. 

The top panel of Fig. 14 also shows that for a cluster 
with a high remnant retention fraction, the impact of the 
disruption time on the MF evolution is similar to that of 
clusters with low retention fractions (see Fig. 9). However, 
the influence of the disruption time becomes smaller when 
more remnants are retained. This explains why Baumgardt 
& Makino (2003) only found a very weak dependence of 
the evolution of the MF on the disruption time (also see 
Fig. 6), since they neglected remnant kick velocities and 
retained all remnants in their simulations. 

Analogous to Fig. 10 in Sect. 5.1, the evolution of the 
MF slope in different mass ranges is shown in Fig. 15 for 
the clusters with t = 1 Myr from Fig. 14. 13 The kick 



13 For the clusters with relatively long disruption times that 
are considered in this section, the variable stellar mass range 
that was introduced in Sect. 5.1 to trace the relation between 
MF slope and M/L ratio gives an evolution of the slope that is 
comparable that for the fixed mass ranges. It is omitted from 
the figures in this section to improve their clarity. 
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Fig. 16. Influence of the black hole retention fraction on the 
M/Lv ratio evolution for an initial cluster mass M; = 10 6 M Q . 
Shown is the relative M/Ly ratio decrease with respect to the 
value expected for stellar evolution (M/Lv) s tev versus the re- 
maining cluster mass fraction. The solid, dashed and dotted 
lines represent <7 kicki bh = {40,80,200} km s~\ corresponding 
to / ret ,bh = {0.219, 0.041, 0.003} for a 10 6 M cluster. 

velocity dispersion has an effect on the MF that is more 
uniform than the consequences of variations in the dis- 
ruption time, leading to very similar slope evolutions in 
the two different stellar mass ranges. Independent of the 
mass range, an increase in remnant retention fraction is re- 
flected by an increase of a. The model that is displayed for 
ffkick,bh = 40 km s _1 , t = 1 Myr, and M ; = 10 6 M Q (the 
upper dashed and solid lines in Fig. 15) marks the transi- 
tion between an increase or decrease of the MF slope by 
dynamical evolution. For an initial / ro t,bh < 0.25, low-mass 
stars are preferentially ejected during most of the cluster 
lifetime, while for / ro t,bh > 0.25 mainly the massive stars 
escape. For shorter disruption times, the transition is lo- 
cated at a smaller black hole retention fraction. 

Because the black hole retention fraction affects the 
overall slope of the stellar MF, the changes in a are matched 
by corresponding changes in the M/L ratio. In Fig. 16, the 
relative M/Ly ratio change due to dynamical evolution is 
shown for same clusters as in Fig. 15. Contrary to the clus- 
ters with low remnant retention fractions in Sect. 5.1, the 
M/L ratio of the clusters in Fig. 16 does not monotonously 
decrease. Close to total disruption, the massive remnants 
are the last bodies to be ejected. During that short phase of 
cluster evolution, the M/L ratio is increased by dynamical 
evolution and exceeds the value it would have due to stellar 
evolution alone. 

The behaviour of M/L ratio for different black hole kick 
velocity dispersions has interesting implications for the re- 
lation between stellar MF slope and M/L ratio, which is 
shown in Fig. 17. In combination with Fig. 12 (note the 
different axes) , it shows possible evolutionary tracks of star 
clusters in this plane, indicating that nearly every loca- 
tion may be reached. However, when limiting ourselves to 
long-lived clusters, Fig. 17 illustrates that these clusters will 
follow a trend of decreasing slope with decreasing M/ L ra- 
tio, albeit with excursions to high M/L ratios and slightly 
higher a close to their total disruption. This explains the 
trend that was found by Kruijssen & Mieske (2009), who 



Fig. 17. Influence of the black hole retention fraction on the 
combined evolution of the MF slope a and the M / Lv ratio for an 
initial cluster mass Mi = 10 6 Mq . Shown is a versus the relative 
M/Lv ratio decrease due to dynamical evolution. All clusters 
start at the vertical line (M/Ly )/(M/Ly ) s tev = 1- Solid lines 
denote the slope in the mass range 0.1 < m/M < 0.5 and the 
dashed lines designate the mass range 0.3 < m/M Q < 0.8, with 
from right to left a kick>b h = {40, 80, 200} km s _1 , corresponding 
to /rct.bh = {0.219,0.041,0.003} for a 10 6 M© cluster. 



considered the relation between the observed MF slopes 
and M/L ratios of Galactic globular clusters. 

The colour change due to dynamical evolution is only 
very small for clusters with tjj°* al > 1.5 Gyr (see Sect. 5.1). 
Because clusters in which remnants are retained are mas- 
sive, their lifetimes are correspondingly long. As a result, 
the colour evolution is largely unaffected for the clusters 
in which the remnant retention fraction could play a role 
(A(V — I) < 0.03 mag). The colour change is even smaller 
if more massive remnants are retained, because then the 
stellar MF more closely resembles its initial form (see the 
upper panel of Fig. 14). Long-lived clusters generally ap- 
pear ~ 0.005 mag bluer in V — I due to dynamical evolu- 
tion during the last ~ 3 — 20% of their lifetimes and reach a 
similar reddening upon their total disruption, which is well 
within observational errors. The colours of old clusters are 
thus only marginally affected by dynamical evolution. 

The evolution of the total remnant mass fraction is 
shown in Fig. 18 for different black hole kick velocity dis- 
persions. The seemingly counterintuitive result is that the 
fraction of the cluster mass that is constituted by remnants 
is smaller when more black holes are retained. As shown in 
Fig. 14, the retention of black holes suppresses the deple- 
tion of the low-mass end of the MF due to the 'sweet spot' 
ejection (see Sect. 3) of massive (~ 1 Mq) stars by the black 
holes. After ~ 1 Gyr, white dwarfs and neutron stars have 
masses that are similar to those of the massive stars, im- 
plying that their ejection rate is also increased when more 
black holes arc retained. Because the total mass constituted 
by white dwarfs and neutron stars is larger than the com- 
bined mass of all black holes, the fraction of the total cluster 
mass that is constituted by remnants decreases if these low- 
mass remnants are ejected by the more massive black holes. 
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Fig. 18. Influence of the black hole retention fraction on the 
total remnant mass fraction. Shown is the ratio of the total 
mass in stellar remnants M SI to the cluster mass M versus the 
remaining cluster mass fraction. The solid, dashed and dotted 
lines represent <7kick,bh = {40, 80, 200} km s _1 , corresponding to 
/ret.bh = {0.219, 0.041, 0.003} for a 10 6 M Q cluster. 

6. Discussion and applications 

The results of this paper show that the stellar MFs in star 
clusters differ strongly from their initial forms due to dy- 
namical cluster evolution. The specific kinds of these dif- 
ferences depend on the properties of the star clusters and 
their tidal environment, most importantly on the disruption 
time, remnant retention fraction, and IMF. 14 

A physical model for the evolution of the stellar MF is 
presented in which two-body relaxation leads to a stellar 
mass dependence of the ejection rate. For any particular 
stellar mass, the ejection rate is determined by the typ- 
ical proximity of that mass to the escape energy and by 
the timescale on which the two-body relaxation with the 
other stars takes place. Combined with a prescription for 
stellar evolution, stellar remnant production, and remnant 
retention using kick velocity dispersions, this provides a 
description for the total evolution of the MF. This descrip- 
tion is independent of the adopted total mass evolution. The 
model shows that the slope of the mass function is a possi- 
ble indicator for the mass fraction that has been lost due to 
dissolution, provided that the IMF does not vary and the 
remnant retention fraction has been fairly similar for young 
globular clusters. 15 

For the exact same initial conditions, the model shows 
excellent agreement with TV-body simulations of the evolv- 
ing MF by Baumgardt & Makino (2003). However, an im- 
portant advantage of the presented model compared to the 
(more accurate) TV-body simulations is its short runtime 
and corresponding flexibility. It can be easily applied to 
compute the evolution of clusters for a large range of initial 
conditions. The results can then be used to identify interest- 
ing cases for more detailed and less simplified calculations 
with TV-body or Monte Carlo models. 

14 Although not specifically shown in this paper (but not sur- 
prisingly) , the differences also depend on the initial-final stellar 
mass relation. 

Any variability of the retention fraction would induce sub- 
stantial scatter, see Sect. 5.2 and Fig. 19. 



The most important simplification of the model is ne- 
glecting the effect of binary encounters on the stellar mass 
dependence of the ejection rate. To incorporate binaries, a 
conclusive census of the binary population in star clusters 
would be required, which is not yet available. Nonetheless, 
it is possible to make a qualitative estimate for the effect 
of binaries. The encounter rate of binaries would typically 
be higher than that of individual stars, because the cross 
section of binaries is larger. This would increase the relative 
escape rate at the stellar mass for which the binary frac- 
tion 16 peaks. This binary fraction is found to increase with 
primary mass (see e.g. Kouwenhoven et al. 2009). Because 
massive stars are removed by stellar evolution, this implies 
that the binary fraction decreases with age, which is in 
agreement with the low binary fraction observed in glob- 
ular clusters (<~ 2%, e.g. Richer et al. 2004). The effect of 
binaries on the evolution of the mass function would thus 
be most notable if the majority of the dynamical mass loss 
occurs at ages < 50 Myr (the typical lifetime of an 8 M 
star), in which case it would somewhat enhance the relative 
escape rate of the most massive stars. 

The model is applied to investigate the influence of the 
disruption time and remnant retention on the evolution of 
the MF and integrated photometric properties of star clus- 
ters. For total disruption times t^°* al < 1 Gyr, the modeled 
relative ejection rate is highest at a certain 'sweet spot' 
mass that is typically 15 — 20% of the mass of the most 
massive objects in the cluster. For longer lifetimes, the evo- 
lution of the MF is dominated by low-mass star depletion, 
unless the retention fraction of massive stellar remnants is 
larger than 0.25. Only in the particular case of such a high 
retention fraction, the M/L ratio is increased by dynami- 
cal evolution when the cluster approaches total disruption. 
In all other scenarios, the M/L ratio decreases because the 
most massive (luminous) stars are kept. 17 When defining 
the slope of the MF in the range 30 — 80% of the maxi- 
mum stellar mass, this gives a clear relation between the 
MF slope and the M/L ratio. For slopes that are defined 
in fixed mass ranges, there is not necessarily a correlation 
between slope and M/L ratio if id°* al < 1 Gyr. In clus- 
ters with a longer total disruption time, both quantities 
are related. Dynamical cluster evolution is found to induce 
some reddening of the integrated cluster colours, amount- 
ing up to 0.1 — 0.2 mag in V — I for total disruption times 
i)£* al < 1.5 Gyr. The fraction of the cluster mass that is con- 
stituted by remnants surprisingly decreases if more black 
holes are retained, because the black holes preferentially 
eject bodies around the masses of white dwarfs and neu- 
tron stars, which contain most of the total remnant mass. 

Contrary to what is suggested by other studies (e.g. 
Baumgardt & Makino 2003; Anders et al. 2009), the evo- 
lution of the MF is not homologous. The reason that these 
studies concluded that its evolution is very similar for all 
clusters (also see Figs. 6 and 7), is that they assumed that 
all remnants were retained. It is illustrated in Fig. 14 that 
the differences between clusters with dissimilar disruption 
times disappear when the retention fraction increases. For 
realistic retention fractions, differences do arise. If two clus- 
ters with different initial masses have the same total dis- 



The fraction of stars residing in binary or multiple systems. 
17 This process differs from a possible variability of the propor- 
tionality between the velocity dispersion and the cluster mass, 
which concerns a much shorter timescale (e.g. Boily et al. 2009). 
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ruption time, their MF evolution will be dissimilar due to 
their different remnant retention fractions and the impact 
of the retained remnants on the dynamical cluster evolu- 
tion. Alternatively, if two clusters have equal initial masses 
but different total disruption times, for instance due to dif- 
ferences in their galactic location or environment, their MF 
evolution will be dissimilar due to the dynamical impact of 
the evolution of the maximum stellar mass. 

The larger variation of MF evolution that is found with 
presented model may also be able to explain observations 
of globular clusters in which the MF cannot be charac- 
terised by a single power law (De Marchi et al. 2000). If 
the evolution of the MF were homologous, these features 
would likely be primordial (Baumgardt & Makino 2003), 
but this is not necessarily the case when using realistic re- 
tention fractions. Most other differences between the results 
presented in Sect. 5 and those from Baumgardt & Makino 
(2003) are also due to their assumption of full remnant 
retention. For example, their M/L ratio evolution shows 
a smaller decrease than in Fig. ff. This is explained in 
Fig. 16, where it is shown that dynamical evolution reduces 
the M/L ratio by a smaller amount if the retention fraction 
is larger. 

Studies on the fractal nature of cluster formation show 
that star clusters are initially substructured (Elmegreen 
2000; Bonnell et al. 2003). Even though this substructure 
is typically erased on a crossing time, it can induce pri- 
mordial mass segregation in star clusters (McMillan et al. 
2007; Allison et al. 2009). The influence of primordial mass 
segregation on the evolution of the MF has recently been 
investigated by Baumgardt et al. (2008) and Vesperini 
et al. (2009). While Baumgardt et al. (2008) do not in- 
clude stellar evolution and concentrate on two-body relax- 
ation, Vesperini et al. (2009) do include stellar evolution. 
They show that for some degrees of primordial mass segre- 
gation, the mass loss by stellar evolution can induce addi- 
tional dynamical mass loss that strongly decreases the to- 
tal disruption time. For clusters that survive for a Hubble 
time, the MF evolution in the case of primordial mass seg- 
regation is very similar to an initially unscgregated clus- 
ter. Vesperini et al. (2009) conclude that the evolution of 
the MF is only affected by primordial mass segregation for 
clusters in which the total disruption time is sufficiently 
decreased by the induced mass loss. In that case, the slope 
of the MF remains much closer to its initial value than 
it would in clusters without primordial mass segregation. 
Their conclusion is consistent with the model presented in 
this paper, because the evolution of the MF is determined 
by the most massive stars at the time when the largest mass 
loss occurs (see Figs. 4 and 9). This induced mass loss en- 
ters the model in terms of the absolute mass loss rate in 
Eq. 3, not in the stellar mass-dependent escape rate per 
unit mass loss rate of Eq. f 3. 

A change in total mass loss rate is not the only con- 
sequence of primordial mass segregation. Baumgardt et al. 
(2008) have shown that low-mass star depletion is enhanced 
for clusters without stellar evolution that are primordially 
mass-segregated. This occurs because energy equipartition 
is reached on a shorter timescale and because of their use 
of a fixed (m max = 1.2 M Q ) maximum stellar mass. As 
a result, there are no massive bodies to increase the ejec- 
tion rate of intermediate mass stars (see Fig. 5), implying 
that only the low- mass stars arc preferentially lost. In the 
present paper, mass segregation is assumed to arise dynam- 
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Fig. 19. MF slope versus remaining lifetime (assuming a glob- 
ular cluster age of 12 Gyr). Diamonds represent the observed 
values from De Marchi et al. (2007), with typical errors as 
shown by the error bar in the lower right corner. The re- 
maining lifetimes are taken from Baumgardt et al. (2008). 
Dotted curves represent the model evolutionary tracks of clus- 
ters with log (Mi/M Q ) = {6,6.25,6.5,6.75,7} from Sect. 5.2 
with {fT k ick,wd,Ckick,n B ,<Tkick,bh} = {4,100,200} km s -1 , corre- 
sponding to {/rct.wd, /ret.ns, /rct,bh} = {0.983, 0.022, 0.003} for a 
10 6 Mq cluster. The solid line connects the present-day loca- 
tions of the modeled clusters in the diagram (crosses), while the 
dashed line represents the same relation for Ukick.bh = 40 km s _1 
(/ret.bh = 0.219 for a 10 6 M cluster). The dash-dotted line 
shows the homologous cluster evolution from Baumgardt & 
Makino (2003). 



ically, but the model could in principle be adapted to cover 
primordial mass segregation by setting C2 = and adjust- 
ing c\ to the initial velocity distribution until it is erased 
by dynamical evolution (see Eq. 22), after which the values 
from Sect. 3 can be used. 18 This does not necessarily yield 
enhanced low-mass star depletion for clusters with a com- 
plete IMF (including masses m > f .2 M ) because of the 
presence of massive stars or remnants. 

The presented model can be applied to the MFs of 
Galactic globular clusters that are observed by De Marchi 
et al. (2007). These MFs are more strongly depleted than is 
found in the A^-body simulations by Baumgardt & Makino 
(2003), which has been attributed to primordial mass segre- 
gation (Baumgardt et al. 2008). However, the observations 
can also very accurately be explained with the realistic rem- 
nant retention fractions that are used in the present paper. 
This is shown in Fig. 19, where the observed MF slopes 
and remaining lifetimes of the globular clusters from De 
Marchi et al. (2007) are compared with the globular cluster- 
like models from Sect. 5.2 (to = 1 Myr). The models are 
in much better agreement with the data than the TV-body 
runs with complete remnant retention from Baumgardt & 
Makino (2003). Deviations to other values of a can occur 
due to variations in disruption time and remnant retention 



18 As explained in Sect 3, ci represents the ratio of the mean 
speed squared to the central escape velocity squared that de- 
pends on the degree of mass segregation (and thus on the IMF). 
On the other hand, C2 is a proportionality constant in the ex- 
pression for the onset of the stellar mass-dependent ejection of 
stars, which depends on the concentration or King parameter. 
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fractions, as is also shown in Fig. 19. For example, a varia- 
tion of the remnant kick velocity with metallicity in combi- 
nation with the known variation of the disruption time (see, 
e.g. Kruijssen & Mieske 2009; Kruijssen & Portegies Zwart 
2009) should be sufficient to cover the observed scatter. 

The above line of reasoning provides an explanation for 
the the depleted MFs in Fig. 19 that is consistent with the 
simulations by Vesperini et al. (2009), who showed that 
the effects of primordial mass segregation are in fact sup- 
pressed in long-lived clusters due to the expansion caused 
by stellar evolution. This increases the relaxation time and 
yields an evolution of the MF that is very similar to the 
initially unscgregated scenario, indicating that primordial 
mass segregation is not a likely explanation for strongly 
depleted MFs. Observations of the remnant composition of 
these globular clusters could reveal a definitive answer as 
to whether the depleted MFs are explained by primordial 
mass segregation or by dynamical evolution with a realistic 
remnant retention fraction. 

Dynamical cluster evolution does not appear to have a 
large effect on the colours of old (globular) clusters. The 
only way in which the colours could be affected beyond 
typical observational errors, is if globular clusters have lost 
substantial fractions of their masses during the first ~ Gyr 
after their formation. In that case, the dynamical evolution 
of the stellar MF in globular clusters may have implications 
for studies of colour bimodality (e.g. Larsen et al. 2001) or 
the blue tilt (e.g. Harris et al. 2006). It could then also 
possibly explain the trend of increasing V — K colour with 
decreasing M/Ly ratio found by Strader et al. (2009) for 
globular clusters in M31, because quickly dissolving clus- 
ters generally become redder and have reduced M/L ratios. 
More research is needed to determine the role of the chang- 
ing MF in the above properties of globular cluster systems. 

It can be concluded that the evolution of the stellar MF 
in star clusters is not as similar for all clusters as previously 
thought. Its precise evolution is determined by cluster char- 
acteristics like the disruption time, the remnant retention 
fraction, initial- final stellar mass relation, and the IMF. In 
order to decipher the evolution of observed star clusters, 
it is essential to record these characteristics and to relate 
them to possible scenarios for the internal evolution of clus- 
ters. That way, observables like the slope of the MF, the 
M/L ratio, the broadband colours, and the mass fraction 
in remnants can be better understood. 
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